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Abstract 

We report on the general purpose numerical program Mercutio, which can be 
used to calculate any infrared safe four-jet quantity in electron-positron annihila- 
tion at next-to-leading order. The program is based on the dipole formalism and 
uses a remapping of phase-space in order to improve the efficiency of the Monte 
Carlo integration. Numerical results are given for the four-jet fraction and the D- 
parameter. These results are compared with already existing ones in the literature 
and serve as a cross-check. The program can also be used to investigate the internal 
structure of three-jet events at NLO. We give results for previously uncalculated 
observables: the jet broadening variable and the softest-jet explanarity. 
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1 Introduction 



Electron-positron annihilation provides a clean environment for studying jet physics, as 
strong interaction effects occur only in the final state. We thus avoid the complications 
of parton distribution functions and their associated uncertainties. Perturbative QCD 
corrections to jet variables can nonetheless be sizeable. The most recent addition to the 
list of known one-loop amplitudes are those for e + e~ — > 4 partons, which have been cal- 
culated by two groups independently. Our group obtained analytic expressions for the 
amplitudes [[J by calculating helicity amplitudes using the cut technique. The Durham 
group performed a conventional calculation and provided computer code . The two sets 
of expressions have been shown to agree numerically. Here we report on the implementa- 
tion of the one-loop amplitudes into a numerical program Mercutio, which can be used 
to calculate any infrared safe four-, three- or two-jet quantity at next-to-leading order 
(and five-jet quantities at leading order). It is similar to Menlo Parc ||, Debrecen gj] 
and EERAD2 [Q. The main difference among the programs is the technique used for the 
cancellation of real and virtual singularities. Menlo Parc uses the subtraction method 
of ref. 0, Debrecen the dipole method 0, while EERAD2 uses a hybrid of the phase- 
space slicing || H and subtraction methods. Mercutio also uses the dipole formalism, 
but with a remapping of phase space. A naive implementation of the dipole formalism 
will give large statistical errors when performing a Monte Carlo integration over the real 
corrections with dipole factors subtracted. In order to improve the efficiency of the Monte 
Carlo integration we remap the phase space in order to make the integrand more flat. 

With the resulting numerical program, we first compute results for several observables 
studied previously in the literature. We present results for the four-jet fraction and the 
D-parameter, which agree with the results from the other programs. We also present re- 
sults for quantities describing the internal structure of three-jet events. We give numerical 
results for the jet broadening and the softest-jet explanarity. 

The paper is organized as follows: in section 2 we explain the implementation of the 
dipole formalism and the remapping of phase space. Section 3 deals with technical de- 
tails concerning various dimensional regularization schemes. Section 4 is devoted to phe- 
nomenology, and contains our numerical results. The conclusions are given in section 
5. 



2 Cancellation of Infrared Divergences 

The NLO cross section receives contributions both from virtual corrections and from real 
emission. For e + e~ — > 4 jets the virtual part consists of the interference terms between the 
tree amplitudes and the one-loop amplitudes for e + e~ — > qggq and e + e~ — > qqq'q'. The 
real emission part consists of the squared tree-level amplitudes with one additional parton, 



namely e + e~ — > qgggq and e + e~ — > qqq'q' g, which were computed long ago [TO]. Both the 
virtual and the real emission contributions contain infrared divergences, arising when one 
particle becomes soft or two particles collinear. Only the sum of the virtual corrections 
and the real emission part is finite. Because the virtual part is integrated over n-particle 
phase-space, whereas the real emission part is integrated over the (n + l)-particle phase 
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space, it is not possible to integrate the divergent contributions numerically and to cancel 
the divergences after the (numerical) integration. The cancellation of infrared singular- 
ities has to be performed before any numerical integration is done. There are basically 
two approaches to overcome this obstacle: phase-space slicing || and the subtraction 
method. (Within the subtraction method there are two variants: the subtraction method 
by Frixione, Kunszt and Signer || and the dipole formalism of Catani and Seymour [0.) 
As mentioned above, Mercutio uses the dipole formalism in order to handle infrared 
divergences. 

Throughout this paper the conventional normalization of the color matrices is used: 
TrT a T b = i S ab_ 

2.1 The dipole formalism 

The dipole formalism [7[] is based on the subtraction method. The NLO cross section is 
written as 

a NLO = J da R + Jda v 

n+1 n 

= J (da R - da A ) + J (da v + J daA , (1) 

n+1 n \ 1 / 

where in the last line an approximation term da A has been added and subtracted. The 
approximation da A has to be a proper approximation of da R such as to have the same 
pointwise singular behaviour in D dimensions as da R itself. Further da A has to be analyti- 
cally integrable in D dimensions over the one-parton subspace leading to soft and collinear 
divergences. According to Catani and Seymour the approximation term is written as the 
n-parton Born term da B times the dipole factors, 

da A = £ da B ®dV d]pole . (2) 
dipoles 

The integration over the (n + l)-parton phase space yields 



da A = Yl I d(yB ® I rfV dipole = / [ d(yB ® J J > ( 3 ) 

n+1 dipoles n 1 n 

where the universal factor I is defined by 

1 = E / ^dipoles- ( 4 ) 
dipoles l 

For the explicit formulae we refer to the paper of Catani and Seymour J7J. Here, we 
discuss the implementation of the dipole formalism when using helicity amplitudes. Let 
denote the n-parton amplitude, where the polarization vector e u of the emitter gluon 
has been amputated. Using 

(4)*£ + + (£ M )*^ = _^ + ^VL_IiVL ( 5 ) 

Pij ■ q 
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where q is an arbitrary reference momentum, we obtain 



{A*y{-g^,)A v = A* + A + + A*_A_ (6) 

where A± denotes the helicity amplitude where the emitter gluon has positive or negative 
helicity. The dependence on the reference momentum q drops out, since (pij) A 1 = 

(gauge invariance), q^ = (property of the polarization vectors) and (pij) 2 = (the 
gluon is on mass-shell). For the spin correlation we obtain 

(A,)* (z iP ? - zffl) (z iP » - z»p») A v = \EA + + E*A^\ 2 , (7) 

where 

{q+\zii>i-z j i>j\p ij +) 
= /or - i W 

vernal 

Using the fact that the spin correlation tensor is orthogonal to pij one shows again that 
the dependence on the reference momentum drops out. As reference momentum one may 
choose q = Pj, in that case E reduces to 

E = ZiipjPiKpipjj) ^ 
\~-i>ji'i.i. 



2.2 Improving the Dipole Formalism 

Within the dipole formalism we have to evaluate the terms 

J (da v + da B ®l) and J da R - da A . (10) 

n n+1 

The first term is most efficiently integrated by splitting it into leading-color and subleading- 
color pieces. The leading-color piece gives the numerically dominant contribution. The 
leading-color one-loop amplitudes are relatively simple and can therefore be integrated 
without further problems. The subleading-color one-loop amplitudes are more compli- 
cated and therefore need more computer time. On the other hand, they will yield contri- 
butions which are numerically smaller. We can thus evaluate them using fewer integrand 
evaluations. Although the second term da R — da A usually gives only a modest numerical 
contribution, a naive Monte Carlo integration will give large statistical errors. Due to 
the large number of dipole factors involved, it is also very computationally intensive. We 
first deal with the fact that we have to reproduce logarithms numerically. This is solved 
by a remapping of phase-space. In a second stage we use color decomposition in order to 
reduce the number of dipole factors that have to be evaluated. 

A simplified model for the term da R — da A would be 

v 7 
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where /(O) = g(0) is assumed. The integral can be rewritten as 

2/min £ { \ { \ ^ 

F = J dx nX) ' x 9[X) + J dy(f(e*)-g(#)), (12) 

In j/ min 

where y min is an artificial parameter separating a numerically dangerous region from a 
stable region. Using the Taylor expansion for f(x) — g(x), one sees that the first term 
gives a contribution of order 0(y ni \ n ). In the second term the 1/x behaviour has been 
absorbed into the integral measure by a change of variables y — \nx, and the integrand 
tends to be more flat. This should reduce the statistical error in a Monte Carlo integration. 

The precise remapping of phase space for the term da R — da A is done as follows: 
Consider a set of products of invariants, including all the invariants in which the (unsub- 
tracted) matrix elements has poles (and which may give rise to logarithms therefore). The 
relevant set of products of invariants for the two-quark, three-gluon final state consists of 
the pairs 

<Sg = { s ql s 12, S q iS 13 , S q2 S 2 i, S q2 S 23 , S q3 S 3 i, S q3 S 32 

S23S35, s 32 s 2q , s 13 s 3q , S U S lq , S 12 S 2q , S 2 iS lq 

s 12 s 23> s 13 s 32, S21S13 

SqlSlq, S q2 S 2q , S q2 S 2q } , (13) 

where we have labelled the gluons from 1 to 3, q corresponds to a quark and q to the 
anti-quark. Singularities associated with the last three pairs appear only in the color- 
subleading part. The corresponding set for the four-quark, one-gluon final state consists 
of the pairs 

<Sq = { s qg s gq-> s q'g s gq' •> s qg s gq' ■> s q'g s gq 

SqqSqq'i Sqq'Sq'q'i ^qq^qq') Sqq'Sq'q' 

s qg s gq'-> s qg s gq'} • (14) 

The second line takes care of the collinear singularities when two quarks become collinear. 
The phase space is then partitioned according to 

$ n+i = ]Te(a,6,c)$ n+ i (15) 
s 

where Q(a,b,c) = 1 if s a bSb c is the smallest product in the set S, and Q(a,b,c) = oth- 
erwise. The sum is over all products in the set S. We may use the symmetries under 
permutations in order to reduce the number of summands which need to be evaluated. 
For the two-quark, three-gluon case we have to evaluate the terms with Q(q, 1,2) and 
0(2, 3, q) with weight 6, as well as the terms with 6(1,2,3) and Q(q,3,q) which are 
weighted by a factor of 3. For the four quark, one gluon case we evaluate the terms with 
Q{q, g,q'),Q(q, g,q), ®(q,q,q') and Q(q,q,q'), which are weighted by a factor 2, and the 
terms with Q(q,g,q') and Q(q,g,q') with unit weight. 

Next, introduce a parameter s m i n which separates numerically stable and unstable 
regions. (It is analogous to the s m i n parameter in the phase-space slicing method, but 
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plays a somewhat different role here.) Using this parameter, define two regions as follows. 
Suppose s as s sb is the smallest product in the set. Then the first region is defined as the 
region where s as > s min and s sb > s m[n . In this region we perform a remapping of the 
phase space as explained in detail in the next subsection. 

The second region is the complement of the first: s as < s min or s sb < s min . In this 
region the integration is performed without any phase space remapping. 

In summary, the numerical computation of the term dcr R — da A goes as follows, 

da R -da A = J2 { daR ~ dcjA ) ( s ^ ~ s min)©(s S 6 - s min )6(a, s, b) 
s 

+ 2 { d(jR ~ dcjA ) ( X ~ ( Sas ~ s min )Q(s sb - s min )) 6 (a, s, b) . 

(16) 

Note that this separation into two regions does not involve any approximations and is 
exact whatever value s m i n might take. The aim is of course to choose s m j n so as to reduce 
the statistical errors. It turns out that choosing s m i n small enough makes the contributions 
of the second region negligible. The contribution of the first region has a statistical error 
reduced by roughly a factor of ten compared to the result without any remapping and 
the same number of integrand evaluations. Empirically 

77 = Smin/Q2 = 10- 5 (17) 

Vent 

is a good value, where y cut is the jet-defining parameter. 



2.3 Remapping of Phase Space 

The remapping of phase space we describe here was originally intended for use with the 
phase space slicing method, and was used in prior jet programs of one of the authors. 
The idea is to generate an (n + l)-parton configuration that is 'near' the hard n-parton 
configuration corresponding to the leading-order calculation. We want to do this in a way 
that makes the (n + l)-parton configuration approach the hard configuration in singular 
limits. 

Suppose we have a set of product of invariants S such that s as s s b is the smallest 
product in the set. In this region we remap the phase space as follows. Let k' a , k s 
and k' b be the corresponding momenta such that s as = (k' a + k s ) 2 , s sb = (k' b + k s ) 2 and 
s a b = (k' a + k s + k' b ) 2 . We want to relate this (n + 1) particle configuration to a nearby 
"hard" n-particle configuration with (k a + k b ) 2 = (k' a + k s + k' b ) 2 , where k a and k b are the 
corresponding "hard" momenta. Using the factorization of the phase space, we have 

dK 2 

d$ n+1 = d^ n _ 1 --d^ 3 (K,k' a ,k s ,k' b ). (18) 
The three-particle phase space is given by 

d$ 3 {K,k' a ,k s ,k' b ) = 1 — ds as ds sb dVt' b d<p s 

S2{2ir)°s ab 

ds as ds sb d(fi s d§2(K,ka,k b ) (19) 



4(27r) 3 s a6 
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and therefore 

Aib Aib ds asds sb d(j> s 

d* n+1 = d* n 4(27r)3sab . (20) 

The region of integration for s as and s sb is s as > s mm , s sb > s min (either from the 0- 
functions of phase space slicing method, or else from the 0-functions separating the two 
regions above) and s as + s sb < s ab (Dalitz plot for massless particles). In order to get a 
flat integrand we want to absorb poles in s as and s sb into the measure. This is done by 
changing the variables according to 

' Q . \ Ul / Q . \ U 2 

13 mm \ / ^ mm 



s s b = s ab (21) 

Sab J V S ab ) 

where < Ui,U2 < 1. Note that Ui,u 2 > enforces s as , s sb > s mm . Therefore this 
transformation of variables may only be applied to invariants Sij where the region < 
Sij < s min is cut out. The phase space measure becomes 

d$ n+1 = d$ n 1 SasSsb In 2 e(s as + s sb < s ab )duidu 2 d(j) s . (22) 

4(2ir) 6 s ab \ s ab J 

This suggests the following algorithm for generating a (n + l)-parton configuration: 

• Take a "hard" n-parton configuration and pick out two momenta k a and k b . Use 
three uniformly distributed ('random') numbers Ui,u 2 ,u 3 and set 

2 

Ml 



Sab 


= (K + fa) 




1 Sxam 


Sas 


= s ab 




V s ab 






Ssb 


= s ab 




V s ab 


<Ps 


= 2nu 3 . 



U2 

1 



(23) 



If (s a s + s sb ) > s ab , reject the event. 



• If not, solve for k' a , k' b and k s . If s as < s sb we want to have k' b — > k b as s as — > 0. 
Define 

jj, S ab S sb ^ S as "T S sb ^ S a b S a s (24) 

2\/ S a b 2^/ S a b %\/ S a b 

9 ab = arccos ( 1 — — — — — ) , 6 sb = arccos (l ) . (25) 

V 2E a E b J V 2E s E b ) v ; 

It is convenient to work in a coordinate system which is obtained by a Lorentz 
transformation to the center of mass of k a + fa and a rotation such that k' b is along 
the positive z-axis. In that coordinate system 

p' a = E a (l, sin 8 ab cos(4> s + ir), sin 8 ab sin(<f) s + ir), cos 6 ab ), 
p s = E s (l,8in0 sb (x>s</> s ,sm6 sb 8m<f> s ,CQs9 sb ), 

P b = E b (l, 0,0,1). (26) 
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The momenta p' a , p s and p' b are related to the momenta k' a , k s and k' b by a sequence 
of Lorentz transformations back to the original frame 



K = ^boostK y {<t>)Kz{0)p' a 



(27) 



and analogous for the other two momenta. The explicit formulae for the Lorentz 
transformations are obtained as follows : 



Let K = yj (k a + k b ) 2 and denote by Pb the coordinates of the hard momentum 
kb in the center of mass system of k a + kb- Pb is given by 



Pb 



""'6 7? i Kb T" 



h • Q 



k° 



K K J " u ' \ K(Q° + K) K 

with Q = k a + kb- The angles are then given by 

6 = arccos 1 - , , 
V WbPb J 

<p = arctan ( — ] 
\PbJ 

The explicit form of the rotations is 

K z {0) = 

Ky{<t>) = 

The boost k' = AboostQ is given by 
k' = 



Q 



f 


1 








\ 







cos 6? 


sin 


e 










1 




V 





— sin (9 


cos 


o J 


/ 


l 








°\ 







COS(/> 


— sin0 










sin0 


COS(/> 





V 











1 / 



q-Q 



+ 



q 



Q° o g ■ Q - 

K q K ,q \K{Q° + K) ' K 



with Q = k a + k b and K = \J [k a + h) 2 . 

• If s as > s s b, exchange a and b in the formulae above. 

• The "soft" event has then the weight 



n+1 — 



2 (2tt) 3 S ab V »„b 

where W n is the weight of the original "hard" event. 



Q 



(28) 



(29) 



(30) 



(31) 



(32) 
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2.4 Color decomposition and soft or collinear limits 

We now turn our attention to the color decomposition of the five-parton tree amplitudes. 
Singling out the numerical dominant contribution, e.g. the leading color part, allows us 
to evaluate this part with fewer dipole factors as subtraction terms. This approach is 
similar to the one followed within the phase-space slicing method. In fact, we have also 
implemented phase-space slicing for the leading color contributions and checked that both 
methods, the dipole formalism and phase-space slicing, give the same numerical results. 
As a side-result we will give the relevant formulae for the contributions from unresolved 
phase-space needed for the phase-space slicing method. 

The color decomposition of the tree-level amplitude for e + e~ — > qgig 2 g 3 q is 

A 5 (q, 01, <fe, <fe, q) = {T l T 2 T\ q A^ {q, 1, 2, 3, q) + (T l T*T 2 ) q ,A {2) (q, 1, 3, 2, q) 

+ (T 2 T 3 T 1 ) gg -4 3) (g, 2, 3, 1, q) + {T 2 T l T\ q Af\q, 2, 1, 3, q) 
+ (T 3 T 1 T 2 ) gg -4 5) (g,3,l,2,g) + (T^T 2 T 1 ) q qAf^ (q, 3, 2, 1, g).(33) 

The index i of A§ labels the permutation of the gluons. Explicit formulas for the partial 
amplitudes are given in ref. ||10|j . 

The tree-level amplitude for e + e~ — > gqqq'q' can be written as 

Ab(9, q, q, q', q') = ^ U T 32 D 1 + ^S 32 T U D 2 - ^S 12 T U D 3 - ^5 M T 12 D 4 (34) 

where the factor 1/2 is due to the conventional normalization of T a and 
D 1 = C o(l)S 1 (0;l,2;3,4) + c (3)S 2 (0;3,4;l,2) 

+S f iav^- (cb(l)S 3 (0; 1, 4; 3, 2) + c (3)£ 4 (0; 3, 2; 1, 4)) , 
D 2 = c (l)B 2 (0; 1,2; 3, 4) + c (3)5 1 (0; 3, 4; 1,2) 

+5 flavor (c (l)S 4 (0; 1, 4; 3, 2) + c (3)£ 3 (0; 3, 2; 1, 4)) , 

Dz = 00(1)^3(0; 1,2; 3, 4) + 00(3)^4(0; 3, 4; 1,2) 
+5 flav (co(l)5i(0; 1, 4; 3, 2) + c (3)S 2 (0; 3, 2; 1, 4)) , 

D A = co(l)^S4(0;l,2;3,4) + c(3)-^S 3 (0;3,4;l,2) 

+S flav (co(l)B 2 (0; 1, 4; 3, 2) + Co(3)Si(0; 3, 2; 1,4)). (35) 

Explicit formula? for the partial amplitudes Bi are again given in ref. JTU] . In this equation, 
c (j) denotes a factor from the electro- weak coupling and depends explicitly on the flavor 
of the quark qj, 



c (j) = -Q 9j +v e v^V z (s 56 ) 
s 

s - M 2 z ~+iY z M z 



-Pz(s) = : m , ^ u ■ (36) 



The electron - positron pair can either annihilate into a photon or a Z-boson. The first 
term in the expression for cq corresponds to an intermediate photon, whereas the last 
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term corresponds to a Z-boson. The left- and right-handed couplings of the Z-boson to 
fermions are given by 



Is — Q sin 9 W —Q sin 9 



Vl = —a T~ ' Ui? = a ( 3T ) 

sin tV cos &w cos tV 

where Q and I3 are the charge and the third component of the weak isospin of the fermion. 
For an electron and up- or down-type quarks we have: 





Q 


h 




( % ) 


( 1/2 


d J 




I -1/2 




-1 


-1/2 



Next, consider the soft and collinear limit of any of these partial amplitudes, e.g. A® or 
Bi. In the soft-gluon limit, a partial amplitude factorizes, 

A n+1 ^g-Eik x ab -A n , (38) 

where the eikonal factors are 

Squaring the amplitude we obtain terms like 

V A (2) * A (1) + 4 (1) *4 (2) 

soft gluon helicity A 
= g 2 (Eik( 2 )-Eik( 1 ) + + Eik( 2 ) + Eik«-) (aW*AW + A^*A^) . (40) 

There are three different cases, depending how the soft gluon is inserted into the partial 
amplitude A^li and A^ +l . 



1. Two common legs, e.g. A^ ] +l = A^ ] +l (..., a, s, b, ...) and A {2) +1 = A {2) +1 (..., a, s, b, ...), 

V- 4(2)*, (1) ,4(1)* 4(2) _ 2 o 9 S «b f 4(2)* 4(1) I 4(D* 4(2)^ 

(41) 

2. One common leg, e.g. A { ^ +1 = A^l^..., a, s, b, ...) and A {2) +1 = aJ 2 ]^..., a, s, c, ...), 

V A (2) * A (1) + 4 {1) *4 (2) 
A 

= ^-2 



, . , (a6) [co] H \ ( 4(2)* ,(i) , 

(as)(sb) [cs][sa} [as][sb] (cs)(sa) J V " » ^ " » >/ 



= ^ . 9 + _s«^\ ^(2)^(1) + A (i) M ( 2) ) (49) 
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3. No common legs, e.g. A^ ] +l = A^j^..., a, s, b, ...) and A^it = d, s, c, ...) 



(2) _ A (2) 



(2)* ,,(1) L .(1)* 4 (2) 

, , w (a6) [erf] , [06] (cd) \ , (2) (1) (1) (2) x 

4(2)* 4(1) 1 4(1)* 4(2)' 



^.(2W1) ,,(1)* 4 (2) 
A 

^ 1 (as)(sb) [cs][sd] [as][s6] (cs){sd) i 



9-2 



Sad ^ Sbd Sbc 



Sas^sc Sas^sd ^bs^sd ^bs^s 



In the collinear limit tree-level partial amplitudes factorize like 
A n+1 ^g £ Split_ A (^,^)A n (...,P A ,...) 

A=+/- 



(43) 



(44) 



where P = p a + pb, p a = zP and pt = (1 — z)P . A, A a and \ denote the corresponding 
helicities. The splitting functions are fi~lf : 

= 



Split q +(g + ,g + ) = Split (g ,g' 



Sp\it g+ (g + ,g-) = V2- 
$p\it g+ (g-,g + ) = V2 
Sph\+(y~,#~) = -V2 



(1-*)* 
Vz(ab) 

3 



1 



y/z{l-z)[ab] 

SpliV(g + ,s + ) = v / 2-=i— - 
V 1 — z(ao) 

Split,- = -v / 2 



VI — i[a&] 



SpliV(^ + ,g + ) = v / 2^T 
y y/z{ab) 

Spht-(^,g + ) = -v / 2 



1 - 2 



Split 9+ (g + ,g-) = v / 2 
Split 9+ (g~,g + ) = 



yfz[ab] 
1 - z 



Split fl _(^~,(/" 
Sph\-(# + ,#" 

Sph\-(# + ,3~ 

Sph\+(<T,^ 
Sph\+(<r,# H 
SpliV(^",g" 
Split q+(g + ,q~ 



(ab) 
V2 Z 



(ab) 



Split (g ,q 



Split „-{q + ,q 



-V2 
-V2 
V2 



(1-*)* 

3 

Z2 



1 



y/z{l-z)(ab) 

yl — z[a&] 

- * 
VI — z(ab) 



V2 



\fz{ab) 
1-z 



[ab] 



V2- 



ab] 



(45) 



We have checked numerically that all five-parton tree amplitudes in the program have the 
correct collinear limits. 



2.5 Real Emission with Two Quarks and Three Gluons 

The matrix element squared of the tree-level amplitude for e + e~ — > qgxgzgzq can be 
written as 



1^0?, Si, 32, 33, g) I 2 = Ac A 



(46) 
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where 



A* = 



and the color matrix is given by 



1 4\ q 



\5 

4 



A 
A 
V A 



\q,l,2,3,q) \ 
\q,l,3,2,q) 
\q,2,3,l,q) 
\q,2,l,3,q) 

\q,3,l,2,q) 
\q,3,2,l,q) ) 



(47) 



C 



/cj 3) 


J3) 
c 2 


r (3) 


r (3) 
c 2 


r (3) 
c 3 


r (3) 


\ 


c (3) 

c 2 


c (3) 




c 3 


c 2 


c 3 




r (3) 


r (3) 




c 2 


r (3) 
c 3 


c 2 




c (3) 

c 2 


c {3) 

c 3 


c 2 






c 3 






r (3) 
c 2 


c 3 


J3) 
c 4 


r (3) 


c 2 




U 3) 


r (3) 
c 3 


J3) 
c 2 


c 3 


c 2 


r (3) 


/ 



(48 



where the color factors are 



c (3) - 



c (3) - 

c 2 — 



c (3) - 

Co — 



r (3) 
C 4 



IN 3 (iV 2_ !)3 

27 iV 2 ' 
1\ 3 (N 2 -!) 2 
2) W 2 ' 

1\3 N 2 _ ! 

27 N 2 ' 

1\ 3 (iV 2 - l)(iV 2 + l) 

2) iV 2 



o_(3) _(3) 
zc 3 c 2 . 



(49) 



If one gluon becomes soft or collinear, the matrix element squared reduces to the one with 
one gluon less 



\A 5 



t ( c? ] R dtag (q, l,2,q) c^>R off (q, 1, 2, q) 



A 



c ( 2 ] Roff{q, 1, 2, g) c [ i'R d iag{q, 2, 1, g) 



(2) 



where the color factors are 



J2) _ ^\ 2 (^-l) 2 



A 4 



(50) 



4 2) 



l\ 2 iV 2 -l 
2/ N 



(51) 



The finite i?-factors, which give the contribution from unresolved phase-space within the 
phase-space slicing method, are given by 

R d iag(q, 1,2, q) = 

2 

= ^{2 Cl (3)R 2 (q,l,2,q) 
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-2cf (2R 2 (q, 1, 2, g) - R^q, 1, g) - i^g, 2, g)) 

+2cf (i? 2 (g, 1, 2, g) - ^(g, 1, g) - i^g, 2, g) + R (q, g)) 

2 

= ^ {-4 3) (i? 2 (g, 1, 2, g) + i? 2 (g, 2, 1, g) - 2R 1 (q, 1, g) - 2R 1 (q, 2, g)) 



+2cf (R 2 (q, 1, 2, g) + i? 2 (g, 2, 1, g) - R^q, 1, g) - i^g, 2, g)) 
-c? (i? 2 (g, 1, 2, g) + i? 2 (g, 2, 1, g) - 2i? (g, q)) 
+2 W} • 



(52) 



We have written the contribution from each color factor separately. The last line corre- 
sponds to the contribution from the A(g, q, q, q', q') amplitudes, when one pair of quarks 
become collinear. We have used the notation 



N C 2C F 



q->qg 



Ro(q, q) 

R 1 (q,l,q) = S(q,l) + S(l,q)-i- 



(53) 



Nr. 



+1 



-Trim "mm 



2C F q ^ qg 
C A , 



9^99 



Sql s lq 



2C F ™ 



St* 



(54) 



i? 2 (g,l,2,g) = S(q,l) + S(l,2) + S(2,q)- — 



2C F q ^ qg 



Sql 



+1 



3 mm "mm 



9^99 
C A T 



Sql Su 



'mm "mm 



9^99 



Sl2 S 2 q 



2C F ™ 



S(a, b) 

I a ->bc( z l, z< l) 



= In 



2 / A* 



+ In 



S2q 
2 / "^nb 



+ 2 In 



In 



Safe 



^ c (^ 2 ) + / a % c (^ 2 ) hi f^y 

\ "min / 



(55) 
(56) 
(57) 



la^bc an d Il^bc are the terms of order 0(e°) and 0(e x ) obtained from the D-dimensional 
splitting functions integrated over unresolved phase space. In the conventional scheme 
they are given by 



9 9 9 



L q^qg 



C F (-ln, 2 -^ + (iln 2 , 2 -I + ^)e + 0(, 2 )), 



W = T^Q + ^ + O^ 2 )). 



(58) 
(59) 
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2.6 Real Emission with Four Quarks and One Gluon 

The matrix element squared of the tree-level amplitude for e + e~ — > gqqq'q' can 



as 



where 



\A 5 (g,q,q,q',q')\ 2 = &CD. 



D 



( A \ 

D 3 
V D A J 



Then the color matrix is given by 



C 



1 


Jl) 

c l 














Jl) 








c l 




c? 




c {1) 

c 2 


c (1) 

c 2 


c (1) 





V 


Jl) 
c 2 


Jl) 
c 2 





cS 1} J 



where 

The color decomposition of the four quark amplitude can be written as 



^(W34-^W 3 2)x(l,2,3,4) 



1 



1 



^14^32 - ^^12^34 ) 8fi av x(l, 4, 3, 2) 



where 

X (1,2,3,4) = A 4 (1,2,3,4) + A 4 (3,4,1,2). 
The matrix element squared we may write as 



I^UI 2 = x 



... JO) JO) \ 
•t c l c 2 } -> 

of cf> ) X 



with 



J0)_l^ c 2 -1 
° 2 "4 JVL 



and 



X(l,2,3,4) 
<S /Jot ,x(l,4,3,2) 
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In the soft gluon limit, the functions behave like 



£i(0; 1,2; 3, 4) 
52(0; 1,2; 3, 4) 
5 3 (0;1,2;3,4) 
5 4 (0;1,2;3,4) 



<?-Eik 32 -A(l,2;3,4), 
^•Eik 14 -A(l,2;3,4), 
5 -Eik 3 4-A(l,2;3,4), 
5 -Eik 12 -A(l,2;3,4). 



(69) 



In the limit where one pair of quarks becomes collinear, the amplitudes factorize as 
^(0; 1,2; 3, 4) - g £ Split P (3, 4)A 4 (1, P, 0, 2), 

A 

5 2 (0;1,2;3,4) - ^8^(3,4)^(1,0,^2), 



B 3 (0;l,2;3,4) 
B 4 (0;l,2;3,4) 



0, 

g ]T Split P (3, 4) (A 4 (l, P, 0, 2) + A 4 (l, 0, P, 2)) . 

A 



One of the functions P>i is redundant, since 

B\ + B 2 — B> 3 — B4 = 0. 
The contribution from unresolved phase space is written as 



cPR diag (l,2,3,A) ^Roff (1,2,3,4) 



4 0) i? o// (l,2,3,4) cr^a 5 (l,4,3,2) 



(o) 



(70) 



(71) 



(72) 



with 



Rdiag(l, 2, 3, 4) — 



r (o) 



i?o(l, 4) + R (2, 3) + — (i? (l, 2) + i? (3, 4)) 



! C 2 



+4^- [i? (l, 4) + R (2, 3) + i? (l, 2) + it> (3, 4) - R (l, 3) - i? (2, 4)] } , 



P o/ /(l,2,3,4) 



' 2^ [Roil, 4) + P (2, 3) + R (l, 2) + R (3, 4)] 



r (o) 



+2cf (l + ^ j [i2o(l, 4) + R (2, 3) + R (l, 2) + i? (3, 4) - R (l, 3) - i? (2, 4)] 



(73) 



2.7 Color Correlation 

The color correlation matrices may be obtained by two ways. The first approach is the 
one given by Catani and Seymour. Within this approach one acts with the color charge 
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operators on the n-parton amplitudes. For example 

(q, 1, 2, q\T q ■ T^q, 1, 2, q) = (T^T}' ql A(q, 1, 2, q)*, T^ q ,A(q, 2, 1, g)* 

rpa OTV (rpl'rparpX rparpV rpl\ ( TiT iq -A(q,l,2,q~) \ 

■l q , q -21r{l 11 -11 1 )-[ T ^A(q,2A,q) ) 
Here the color charge operator for the gluon has been written as 

ijcab = 2Tr (r c T a T b — T a T c T b ) (75) 



(74) 



On the other hand one may start from the color decomposition of the (n + 1) -parton 
matrix element in the form of equations ([46"D and (|60|) . One then considers the soft and 
collinear limits using the partial fraction decompositions (41) - (|43|). This procedure is 



identical to the one followed in the phase space slicing approach. With the help of the 
identity 

Sab $ab $ab 

+ — ? ; v ( 76 ) 



Sas^sb Sas^Sas ^sb) ^sb\^as S s b) 

the color correlation matrix can then be read off from equations (|5C|) and ([72]). This 
approach has the advantage that it makes the connection between each divergent term 
in the (n + 1) -parton matrix element and the corresponding subtraction term transparent. 

The leading order matrix element e + e~ — > qg\g29^q needs 27 dipole factors. There are six 
terms where the quark is the emitter and a gluon the spectator and three terms where 
the role of the emitter and spectator are exchanged. The color correlation matrices are 
invariant under the exchange of emitter and spectator. The color correlation matrix for 
the case T q ■ T\ is given by 

„(3) _ _(3) i / (3) _ (3) 
c 2 c l 2 V 4 c 2 

2 ( C 4 ^ _ C 2 



T l' Tl - I 1 /,(3) f3)\ J3) (3) I ( 77 ) 



If the antiquark q is the spectator we obtain 



r (3) _(3) 



Tq ' T ~ q ~ { -ci 3) -4 3) ) ^ 

Finally the color correlation matrix where both emitter and spectator are gluons is given 
by 

/ _ c ( 3 ) + 2c (3) - c (3) \ 

Tl ' T2 = { 2 3 -c?+2cf -of ) (79) 

All other color correlation matrices can be obtained by a permutation of indices. 

The leading order matrix element e + e~ — ► gqqq'q' needs 12 dipole factors associated 
with the splitting g — > qq and 12 dipole factors associated to the splitting q — > qg or 
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q — ► gg. The relevant color correlation matrices for the first case where already given 
above. The color correlation matrices for the latter case are 



T ■ T- 



( 



JX) JX) 



H. 

N c 



— ~tr C 



(1) 



4 15 r (D 



i 



— ci - 2 



,(!) 



SX 



O c 2 



H. _ C W (l + Jj\ 



c (1) 



r (1) 

" JV„2 Z AT, 



(81) 



T 9 " T <?' 



c (1) 
9 L 2 



JX 

c 2 



1 + 



r (l) 



c (1) 
O c 2 

Z N c 



(82) 



We can now pick out partial amplitudes corresponding to one specific color factor. This 
subset will in general require fewer dipole factors as subtraction terms. In practise it 
is enough to single out the leading color contribution for the matrix element e + e~ — > 
Q919293Q- Using the symmetry under the permutation of the gluons it is sufficient to 
evaluate one partial amplitude, which we may take to be 



Ag\q,l,2,3,q), 



(83) 



weighted by a factor of six. This partial amplitude needs only six dipole factors Vij k, 
where (ij) denotes the emitter and k the spectator: 



K?l,2> n.2,g) V\ 



12,3? ^23,1> ^23,?) ^3g,2- 



The color correlation matrices simplify accordingly, for example 



T q -T x 



.(3) 







(85) 



which requires only the evaluation of one partial amplitude of the four-parton tree-level 
amplitude. We may further restrict the set of products of invariants to 



\SqlSi2, S12S23, S23S3g} 



(86) 



since we can only have poles in s q \, su, S23 or s 35 . Singling out the leading color con- 
tributions for e + e~ — > qgxQigzq together with the remapping of phase space discussed in 
section [2.3| improves the efficiency of the Monte Carlo integration of the term da R — da A 
sufficiently so that most of the computer time is spent on the evaluation of the virtual 
one-loop amplitudes. 



We finally remark that Z. Nagy and Z. Trocsanyi M have chosen a different approach 
in order to improve the efficiency of the dipole formalism: They introduce an additional 
parameter a and constrain the phase space over which the dipole factors are subtracted: 



dipoles 



dV, 



dipole 



®{Vij,k < «) 



i7) 
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where the dimensionless variable is given by 

PiPi 



Uij,k 



PiPj + PjPk + PiPk 



The three methods (remapping of phase space, color decomposition and restriction of the 
phase space for the dipole factors) could be combined, if there is any need for further 
improvement in the efficiency 

3 Regular izat ion Schemes and Splitting Functions 

Theoretical calculations of infrared-safe quantities in QCD should lead to unambiguous 
results, independent of the chosen regularization scheme. Terms sensitive to the precise 
definition of the regularization scheme enter the calculation usually in the virtual part 
through tensor loop integrals, and in the real emission part through the splitting func- 



tions [12 1, which enter the dipole factors (in the dipole formalism) or the contribution 
from unresolved phase space ( within phase space slicing). Within dimensional regular- 
ization the most commonly used schemes are the conventional dimensional regularization 
scheme (CDR) fl3f , where all momenta and all polarization vectors are taken to be in D 
dimensions, the 't Hooft-Veltman scheme (HV) []H| |, where the momenta and the helicities 



of the unobserved particles are D dimensional, whereas the momenta and the helicities of 
the observed particles are 4 dimensional, and the four- dimensional helicity scheme (FDH) 



lq| , where all polarization vectors are kept in four dimensions, as well as the momenta 
of the observed particles. Only the momenta of the unobserved particles are continued to 
D dimensions. 

The procedure adopted in the numerical program is as follows. The one-loop amplitudes, 
calculated in the four- dimensional helicity scheme, were converted to the HV scheme, 
using the transition rules fT6|| . 

tone-loop tone-loop 2 Atree ( N c T^q n q N c \ /on^ 

■"-fdh ~A HV - c T g A ^— -—- + -—) (8J) 

where n q is the number of quarks. The one-loop amplitude for e + e _ — ► qqgg is converted 
into the 't Hooft-Veltman scheme by 

Aone-loop _ tone-loop 1 2 at [ i __J_ \ Atree / Qn \ 

A HV - A FDH ~ C T~g iV C II - I A , (9(J) 

whereas the one-loop amplitude for e + e~ — > qqq'q' is converted by 

Aone-loop _ A one-loop 2iir ( % ^ \ *tree / Q1 \ 

a hv - a fdh ~ °rg iV c I - - -^2" I A . (91) 

The splitting functions entering the contribution from the integrals over the dipole factors 
have then to be taken in the 't Hooft-Veltman scheme as well. We take the formulae for 



the splitting functions in the HV-scheme from Catani, Seymour and Trocsanyi [17]. They 
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are given by 

P™ = 2C A (-^— z + ^ + z(l-z)y 
= 2C F ( I ^ + (l-£)(l-z)), 

= 2^(1-^4^(1-*)) (92) 
A statistical factor of 1/2! for two identical particles is included in the g — > gg case. 

4 Phenomenology 

The numerical program Mercutio developed in this work is a general purpose program 
for calculating any infrared safe four-, three- and two-jet observable to next-to-leading or- 
der in a s (and five-jet quantities at leading order). It is written in the language C++. The 
only approximations which have been made are the neglect of the light quark masses and 
terms which are suppressed by \jm\ or higher powers of the top quark mass. Mercutio 
calculates the quantities in fixed order in a s , no resummation of terms of In 2 y cut or lny cut 
has (yet) been implemented. Therefore our results are reliable only for values of y cut which 
are not too small. The program calculates the jet quantities at the partonic level, and no 
hadronization is done. 

The distribution for an infrared-safe variable O at the center-of-mass energy Q at next- 
to- leading order is given by two coefficients, Bo and Co, which represent the leading and 
next-to-leading order perturbative contributions: 

h°% - ©'m^vhsm ,93) 

We normalize all distributions to the total hadronic cross section a tot at 0(a s ), given by 

= °£5s( 1 + v)- (94) 

We further normalize each variable such that it takes values between and 1. The average 
of the variable is then defined as 

<°> " h\°> 

As our nominal choice of input parameters we use N c = 3 colors and Nf = 5 massless 
quarks. We take the electromagnetic coupling to be a(m z ) = 1/127.9 and the strong 
coupling to be a s (m z ) = 0.118. The numerical values of the Z°-mass and width are 
m z = 91.187 GeV and T z = 2.490 GeV. For the top mass we take m t = 174 GeV and 
for the weak mixing angle sin 2 6y/ = 0.230. We take the center of mass energy to be 
y/Q^ = m z and we set the renormalization scale equal to fi 2 = Q 2 . 
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Table 1: The four-jet fraction as calculated by Mercutio, MENLO PARC, DEBRECEN 
and EERAD2, for different jet algorithms and varying y cut . 



4.1 Four- Jet Fraction 

The four-jet fraction has been calculated by each group which has provided a numerical 
NLO four-jet program and serves as a cross-check. The four-jet fraction is defined as 

Ri = <I±±*L. (96) 
O'tot 



The values obtained for the four-jet fraction for different jet algorithms ( JADE 



DURHAM |L9[ and GENEVA (fJUf ) and varying y cut are given in table 1, together with 



the corresponding values from the programs MENLO PARC by L. Dixon and A. Signer 
§, DEBRECEN by Z. Nagy and Z. Trocsanyi |§ and EERAD2 by J.M. Campbell, M.A. 
Cullen and E.W.N. Glover 0. 
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4.2 The D-Parameter 



Global event shape variables describe the topology of an event. They may be calculated 
without reference to a jet defining algorithm. Like jet algorithms they have to be infrared 
safe. An example of a global event shape variable is the D-parameter |H] which is derived 
from the eigenvalues of the momentum tensor 



J2 PaPa 



13 = (97) 

E \Pa\ 
a 

where the sum runs over all final state particles and p\ is the i-th component of the three- 
momentum p a of particle a in the cm. system. The tensor 9 is normalized to have unit 
trace. In terms of the eigenvalues of the 9 tensor, Ai, A2, A3, with Ai + A2 + A3 = 1, one 
defines 

D = 27AiA 2 A 3 = 27det0* J '. (98) 

The range of values is < D < 1. The D-parameter measures aplanarity, since one needs 
at least four final-state particles to obtain a non-vanishing value. 

The values for the functions Bo and Co are given in table 2. Figure 1 shows the D- 
parameter distribution. We remark that we have chosen to normalize the distribution to 
the total hadronic cross-section at 0(a s ). The average of the shape variable is defined as 

1 

1 f „ da 



For the average we obtain 



oc. \ 2 , . , . „i fat, x :; 



(D) = — D—dD. (99) 

(T t ot i dV 







(D) = (5.82 ± 0.01) • 10 1 + (2.43 ± 0.06) ■ 10 3 . (100) 

\2tcJ \2tcJ 

These numbers agree with the ones given by Nagy and Trocsanyi as published in ref. 0] 
after taking care of different normalizations (cr tot in our case and <J2-jet i n re f- H)- The 
D-parameter distribution has also been calculated by Campbell, Cullen and Glover ||. 

4.3 The Jet Broadening Variable 

With the numerical program for e + e~ — > 4 jets one may also study the internal structure 
of three-jets events. One example is the jet broadening variable defined as 



1 EWI 

Bjet = —Y, JL T^T ( 101 ) 

n jets j ets E \Pa I 

Here is the momentum of particle a transverse to the jet axis of jet J, and the sum over 
a extends over all particles in the jet J. The jet broadening variable Bj et considered here 
shall not be confused with the narrow jet broadening B m i n , which is defined differently 
The later one has been recently calculated to NLO by Campbell, Cullen and Glover p. 
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Table 2: The Born level and next-to-leading order functions Bp and Co for the D- 
parameter. 
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Figure 1: The D-parameter distribution at NLO (diamonds) and LO (crosses). 



The jet broadening variable is calculated for three-jet events defined by the DURHAM 
algorithm and y cu t = 0.1. This choice is motivated by a recent analysis of the Aleph 
collaboration p2|. For the average we obtain 



(B jet ) = (g) (1.79 ± 0.01) • 10 1 + (g) (4.87 ± 0.31) ■ 10 2 . (102) 

The values for the functions Bs jat and CB jet are given in table 3. Figure 2 shows the 
distribution of the jet broadening variable. 

4.4 The Softest-jet Explanarity 

For three jet events we define the softest-jet explanarity as follows: 

E Pa ■ (Pi X P 2 

n \ 



Esjet 



E\Pa\ PlXP 2 



(103) 



where p a is the three-momentum of particle a inside the softest jet and the sum runs over 
all particles in the softest jet. Pi and P 2 denote the three-momenta of the remaining two 
jets in the event. The explanarity measures the degree to which the shape of the softest 
jet lies out of the event plane. The definition is motivated by the fact, that the explanarity 
will be sensitive to the color factors of QCD. 

The softest-jet explanarity is calculated for three-jet events defined by the DURHAM 
algorithm and y cut = 0.1. For the average we obtain 

(E 3jet ) = (J^Y (1.29 ± 0.01) • 10 1 + {^\ (3.33 ± 0.22)- 10 2 . (104) 
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B Bjet 








0.02 


(7.31 ±0.22) 


•10 1 


(-1.71 ±0.66) 


• 10 3 


0.06 


(5.13 ±0.05) 


•10 1 


(9.23 ±2.79) • 


10 2 


0.10 


(4.27 ±0.06) 


•10 1 


(1.27 ±0.20) • 


10 3 


0.14 


(3.75 ±0.06) 


•10 1 


(1.24 ±0.16) • 


10 3 


0.18 


(3.36 ±0.06) 


•10 1 


(1.26 ±0.14) • 


10 3 


0.22 


(2.99 ±0.03) 


•10 1 


(1.28 ±0.14) • 


10 3 


0.26 


(2.74 ±0.03) 


•10 1 


(1.13 ±0.12) • 


10 3 


0.30 


(2.50 ±0.02) 


•10 1 


(1.06 ±0.07) • 


10 3 


0.34 


(2.28 ±0.03) 


•10 1 


(1.00 ±0.07) • 


10 3 


0.38 


(2.08 ±0.02) 


•10 1 


(9.12 ±0.51) • 


10 2 


0.42 


(1.87 ±0.01) 


•10 1 


(7.92 ±0.69) • 


10 2 


0.46 


(1.67 ±0.02) 


•10 1 


(7.59 ± 1.19) • 


10 2 


0.50 


(1.44 ±0.02) 


•10 1 


(6.36 ± 1.43) • 


10 2 


0.54 


(1.20 ±0.01) 


• 10 1 


(5.75 ±0.30) • 


10 2 


0.58 


(9.42 ±0.18) 


• 10° 


(4.30 ±0.33) • 


10 2 


0.62 


(6.46 ±0.12) 


•10° 


(3.33 ±0.65) • 


10 2 


0.66 


(3.76 ±0.07) 


• 10° 


(1.84 ±0.18) • 


10 2 


0.70 


(1.60 ±0.04) 


•10° 


(7.97 ± 1.07) • 


IO 1 


0.74 


(3.39 ±0.28) 


io- 1 


(1.95 ±0.55) • 


IO 1 


0.78 


(1.11 ±0.31) 


io- 2 


(2.79 ± 1.24) • 


10° 


0.82 


0.00 




(8.26 ±9.73) • 


io- 2 


0.86 


0.00 




0.00 




0.90 


0.00 




0.00 




0.94 


0.00 




0.00 




0.98 


0.00 




0.00 





Table 3: The Born level and next-to- leading order functions B Bjet and Cs. et for the jet 
broadening variable. 
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Figure 2: The S jet distribution at NLO (diamonds) and LO (crosses) 



The values for the functions B E:j . et and C^ 3jet are given in table 4. Figure 3 shows the 
distribution of the softest-jet explanarity. 



5 Conclusions 

In this paper we reported on the general purpose numerical program "Mercutio" , which 
can be used to calculate any infrared safe four-jet quantity in electron-positron annihi- 
lation at next-to-leading order. The program is based on the dipole formalism and uses 
a remapping of phase-space in order to improve the efficiency of the Monte Carlo inte- 
gration. We presented numerical results for the four-jet fraction and the Z)-parameter. 
These results agree with the results from other groups. The program can also be used to 
investigate the internal structure of three-jet events at NLO. We also presented results 
for previously uncalculated observables: the jet broadening variable and the softest-jet 
explanarity. 
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Figure 3: The E 3 j et distribution at NLO (diamonds) and LO (crosses). 
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